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Numerical  calculations  were  carried  out  to  examine  the  physics  of  the  operation  of  a 
nanosecond-pulse,  single  dielectric  barrier  discharge  in  a  configuration  with  planar  sym¬ 
metry.  This  simplified  configuration  was  chosen  as  a  vehicle  to  develop  a  physics-based 
nanosecond  discharge  model,  including  realistic  air  plasma  chemistry  and  compressible 
bulk  gas  flow.  First,  a  reduced  plasma  kinetic  model  (15  species  and  42  processes)  was 
developed  by  carrying  out  a  sensitivity  analysis  of  zero-dimensional  plasma  computations 
with  an  extended  chemical  kinetic  model  (46  species  and  395  processes).  Transient,  one¬ 
dimensional  discharge  computations  were  then  carried  out  using  the  reduced  kinetic  model, 
incorporating  a  drift-diffusion  formulation  for  each  species,  a  self-consistent  computation 
of  the  electric  potential  using  the  Poisson  equation,  and  a  mass-averaged  gas  dynamic 
formulation  for  the  bulk  gas  motion.  Discharge  parameters  (temperature,  pressure,  and 
input  waveform)  were  selected  to  be  representative  of  recent  experiments  on  bow  shock 
control  with  a  nanosecond  discharge  in  a  Mach  5  cylinder  flow.  The  computational  results 
qualitatively  reproduce  many  of  the  features  observed  in  the  experiments,  including  the 
rapid  thermalization  of  the  input  electrical  energy  and  the  consequent  formation  of  a  weak 
shock  wave.  At  breakdown,  input  electrical  energy  is  rapidly  transformed  (over  roughly 
1  ns)  into  ionization  products,  dissociation  products,  and  electronically  excited  particles, 
with  subsequent  thermalization  over  a  relatively  longer  time-scale  (roughly  10  ps).  The 
motivation  for  this  work  is  modeling  nanosecond-pulse,  dielectric  barrier  discharges  for 
applications  in  high-speed  flow  control.  The  effectiveness  of  such  devices  as  flow  control 
actuators  depends  crucially  on  the  rapid  thermalization  of  the  input  electrical  energy,  and 
in  particular  on  the  rate  of  quenching  of  excited  electronic  states  of  nitrogen  molecules  and 
oxygen  atoms  and  on  the  rate  of  electron-ion  recombination. 


I.  Introduction 

Interest  in  plasma-based  flow  control  dates  to  the  mid-1950s,  when  magnetohydro  dynamic  reentry  heat 
shields  were  first  investigated.1,2  Activity  in  the  research  area  waned  in  the  1970s,  with  some  work  on 
drag  reduction  using  corona  discharges  appearing  in  the  1980s.3  A  resurgence  in  the  field  took  place  in  the 
1990s,  with  the  introduction  of  dielectric  barrier  discharge  (DBD)  actuators,4  a  revisit  of  reentry  magne- 
tohydrodynamics,5  and  the  disclosure  of  the  AJAX  hypersonic  vehicle  concept.6  In  the  fifteen  years  since 
this  resurgence,  plasma-based  flow  control  techniques  have  been  a  topic  of  ongoing  research,  motivated  by 
the  possibility  of  extremely  rapid  actuation,  a  low-profile  configuration,  and  the  ability  to  operate  in  hostile 
environments . 71 7 

In  the  high-speed  regime,  plasma-based  flow  control  devices  have  suffered  the  drawbacks  of  either  ex¬ 
cessive  weight18  or  insufficient  control  authority.4  Pulsed  discharge  devices,  based  on  either  arc19  or  glow 
discharges,20  seem  to  be  a  promising  way  around  these  difficulties.  Nanosecond-scale  pulsed  glow  discharges 
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are  efficient  generators  of  both  ions  and  electronically  excited  molecules  because  of  the  extremely  high  in¬ 
stantaneous  reduced  electric  field.21 

The  generation  of  shock  waves  by  volumetric  heat  release  in  pulsed  discharges  was  observed  and  explained 
in  the  1970s  in  the  context  of  gas  laser  technology.22,23  Early  computations  by  Aleksandrov  et  al.23  assumed 
that  all  the  power  dissipated  in  the  discharge  immediately  went  into  heating  the  neutral  gas.  Popov,25 
however,  has  emphasized  a  two-stage  heating  mechanism  in  which  product  species  and  electronically  excited 
species  are  generated  by  electron  impact,  and  then  the  stored  chemical  energy  is  converted  to  thermal 
energy  through  quenching  reactions.  Two-dimensional  calculations  have  been  carried  out  recently  by  Unfer 
and  Boeuf,24  assuming  instant  thermalization  of  30%  of  the  dissipated  power  going  into  electronic  excitation. 

In  recent  experiments,17,26  control  of  a  Mach  5  cylinder  flow  was  demonstrated  using  a  pulsed  surface 
dielectric  barrier  discharge.  Schematics  of  the  experiment  are  shown  in  Fig.  1.  Figure  la  shows  a  diagram  of 
the  model  mounted  in  the  wind  tunnel,  and  Figs,  lb-c  illustrate  the  actuator  configuration  on  the  cylinder. 

The  hollow  cylinder  model  was  made  of  fused  quartz,  with  a  6  mm  outside  diameter  and  a  2  mm  thick 
wall.  A  thin  copper  exposed  electrode  (12  mm  x  1.5  mm  x  0.1  mm)  was  affixed  to  the  surface  of  the  cylinder, 
with  a  second  copper  electrode  mounted  inside  (a  3  mm  diameter  tube,  0.35  mm  thick  and  10  mm  long). 
DuPont™  Kapton®  polyimide  film  was  placed  over  the  ends  of  the  exposed  electrode,  leaving  a  10  mm 
span  exposed  to  the  flow. 

A  combination  of  positive  and  negative  polarity  pulses  to  the  two  electrodes  produced  a  potential  dif¬ 
ference  of  about  27  kV,  lasting  on  the  order  of  5  ns  (pulse  full- width  at  half  maximum).  The  effects  of 
the  energy  release  in  the  resulting  discharge  were  captured  using  phase-locked  schlieren  imaging.  Side-view 
schlieren  images  are  shown  in  Fig.  2  for  several  time-delays  after  the  discharge.  A  weak  shock  wave  is  seen 
to  form  near  the  edge  of  the  exposed  electrode,  and  propagate  upstream  in  the  shock-layer  flow  on  a  time 
scale  on  the  order  of  microseconds.  When  it  reaches  the  bow  shock,  the  shock  shape  is  altered,  and  the 
shock  standoff  increases  by  up  to  25%.  Given  the  close  relationship  between  shock  standoff  and  gradients 
at  the  stagnation  point,9  this  system  might  be  able  to  alter  heat  transfer  rates  at  the  nose  of  blunt  body 
flows.  Further,  the  system  is  promising  for  general  high-speed  flow  control  applications,  for  example  control 
of  supersonic  inlet  flows.27 

We  have  begun  to  formulate  a  high-fidelity  physical  model  of  the  energy  transfer  process  in  the  pulsed 
surface  dielectric  barrier  discharge.28  For  simplicity  in  the  present  work,  we  have  focused  on  a  planar 
geometry,  where  experimental  evidence  exists  that  a  nearly  one-dimensional  discharge  can  occur  at  relatively 
low  pressures.29  (In  a  companion  project,  we  are  exploring  the  three-dimensional  fluid  mechanics  of  the 
experiment  using  a  gasdynamics  code  with  a  phenomenological  volumetric  energy  deposition  model.30)  Using 
coupled  modeling  of  the  plasma  and  compressible  flow  in  a  one-dimensional  geometry,  we  plan  to  study  the 
dominant  physical  effects,  including  energy  thermalization  kinetics  and  compression  wave  formation  and 
propagation. 

To  this  end,  a  reduced  plasma  kinetic  model  (15  species  and  42  processes)  was  developed  first  by  carrying 
out  a  sensitivity  analysis  of  a  zero-dimensional  plasma  computation  with  an  extended  chemical  kinetic 
model  (46  species  and  395  processes).  Transient,  one-dimensional  discharge  computations  were  then  carried 
out  using  the  reduced  kinetic  model,  incorporating  a  drift-diffusion  formulation  for  each  species,  a  self- 
consistent  computation  of  the  electric  potential  using  the  Poisson  equation,  and  a  mass-averaged  gas  dynamic 
formulation  for  the  bulk  gas  motion. 

II.  Development  of  Reduced  Kinetic  Model 

To  obtain  the  reduced  kinetic  model,  we  applied  a  sensitivity  analysis  to  a  detailed,  transient,  zero¬ 
dimensional  air  plasma  model  used  in  previous  work.31  The  reduced  model  was  developed  to  identify  the 
dominant  species  and  reactions  affecting  the  energy  balance  and  the  rate  of  thermalization  in  the  discharge, 
and  to  minimize  the  computational  cost  of  the  transient,  one-dimensional  calculations  that  will  be  presented 
in  Sec.  III. 

The  full  air  plasma  model  is  based  on  the  model  developed  by  Kossyi  et  al.32  It  incorporates  a  set 
of  ordinary  differential  equations  for  number  densities  of  neutral  species  N,  N2,  O,  O2,  O3,  NO,  NO2, 
N20,  N03,  charged  species  e",  N+,  N+,  N+,  N+,  0+,  O+,  O+,  NO+,  NO+,  N20+,  N20+,  N2NO+, 
02NO+,  NONO+,  O",  02“,  O3  ,  NO-,  NO^T,  NO3  ,  N20",  and  excited  species  N2(A3F),  N2(B3II),  N2(C3II), 
N2(a/:LF),  02(a1A),  02(b1F),  02(c1F),  N(2D),  N(2P),  0(XD)  produced  in  the  plasma,  as  well  as  an  energy 
equation  for  predicting  the  time  evolution  of  gas  temperature.  This  set  of  equations  is  coupled  with  a  steady, 
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two-term  expansion  of  the  Boltzmann  equation  for  the  electron  energy  distribution  function  (EEDF)  of 
the  plasma  electrons.33  These  calculations  employ  experimental  cross  sections  of  electron  impact  electronic 
excitation,  dissociation,  ionization,  and  dissociative  attachment  processes.34,35  The  rate  coefficients  of  these 
electron  impact  processes,  as  functions  of  the  reduced  electric  field  E/N ,  are  derived  from  the  Boltzmann 
solutions  by  averaging  the  cross  sections  over  the  EEDF.  The  model  also  incorporates  chemical  reactions 
of  ground  state  species  and  excited  electronic  species,  electron-ion  recombination  and  ion-ion  neutralization 
processes,  ion-molecule  reactions,  and  electron  attachment  and  detachment  processes.  The  rate  coefficients 
of  these  processes  are  taken  from  Kossyi  et  al.32  The  air  plasma  processes  and  the  kinetic  rates  used  are 
listed  in  a  previous  paper.31 

In  the  zero-dimensional  modeling  calculations  using  the  full  air  plasma  kinetic  model,  the  time-resolved 
reduced  electric  field  (i.e.  the  pulse  voltage  waveform)  was  one  of  the  inputs  for  the  model.  The  voltage 
waveform  was  approximated  as  a  pulse  with  Gaussian  form,  U(t)  =  C/peak  exp[— (t  —  to)2/r2],  where  I/peak  is 
the  peak  voltage,  to  is  the  moment  when  the  voltage  peaks,  and  r  =  3  ns  is  the  pulse  width  parameter.  The 
pulse  width  parameter  was  chosen  to  fit  the  experimental  voltage  pulse  duration  produced  by  a  nanosecond 
pulse  generator  (FID  GmbH,  FPG  60-100MC4)  with  pulse  full-width  half-maximum  of  5  ns.26 

The  pulse  peak  voltage  was  considered  an  adjustable  parameter,  and  varied  to  produce  discharge  pulse  en¬ 
ergy  loading  of  50-100  me V/ molecule.  At  these  conditions,  if  all  discharge  input  energy  were  thermalized,  the 
resultant  temperature  rise  in  the  discharge  would  be  AT  ~  165-330  K  (thermalization  of  1  meV/molecule 
input  energy  corresponds  to  AT  =  3.32  K).  However,  recent  kinetic  modeling  calculations  and  experi¬ 
ments25,36,37  suggest  that  approximately  30%  of  the  input  energy  is  thermalized  rapidly  after  a  discharge 
pulse  with  a  peak  reduced  electric  field  of  (T/7V)peak  =  200-900  Td  (1  Td  =  10-17  V-cm2),  during  collisional 
quenching  of  electronically  excited  species. 

Calculations  and  experiments  by  Aleksandrov  et  al.38  suggest  that  in  nanosecond  pulse  discharges  at 
atmospheric  pressure  at  very  high  values  of  the  reduced  electric  field,  E/N  ~  1000  Td,  the  rapidly  thermalized 
energy  fraction  increases  up  to  about  50%,  due  to  contributions  of  ion- molecule  reactions,  electron- ion 
recombination,  and  ion-ion  recombination.  The  time  scale  for  the  rapid  energy  thermalization  ranges  from 
a  few  microseconds  at  p  ~  0.01  atm37  to  below  one  microsecond  at  p  ~  1  bar.38  This  effect  would  limit  the 
temperature  rise  due  to  rapid  energy  thermalization  after  the  discharge  pulse  to  AT  ~  50-100  K.  This  is 
consistent  with  the  temperature  rise  measured  in  a  single-pulse  nanosecond  surface  dielectric  barrier  discharge 
(SDBD)  in  dry  air  at  p  ~  30  Torr  using  the  FID  pulse  generator,  AT  =  40  db  30  K26  and  in  a  single-pulse 
nanosecond  SDBD  in  room  air  using  a  custom-designed  nanosecond  pulse  generator,  AT  =  80  db  50  K.39 

In  the  present  paper,  our  primary  objective  is  to  study  the  effect  of  rapid  energy  thermalization  on 
compression  wave  formation  in  the  discharge,  for  a  time  scale  after  the  discharge  pulse  shorter  than  the 
acoustic  time  Tacoustic  =  E/a.  (Here  L  is  the  characteristic  size  of  the  plasma  and  a  is  the  local  speed  of  sound.) 
Strong  compression  waves  generated  by  nanosecond-pulse  discharges  can  be  used  for  high-speed  flow  control. 
Nanosecond-gate,  broadband  ICCD  images  of  nanosecond-pulse  surface  dielectric  barrier  discharges17, 26 
have  shown  that  the  thickness  of  the  near-surface  plasma  layer  at  atmospheric  pressure  is  of  the  order  of 
100  pm,  which  corresponds  to  an  acoustic  time  of  rac0ustic  ^  0-3  ps.  Thus,  nanosecond-pulse  discharge 
energy  thermalization  on  a  shorter  time  scale  at  atmospheric  pressure  would  result  in  strong  compression 
wave  formation,  as  detected  in  experiments.17,20,39,40 

Transient,  zero-dimensional  calculations  have  been  conducted  for  a  single-pulse  discharge  in  dry  air  at 
pressures  ranging  from  30  Torr  to  760  Torr,  initially  at  room  temperature.  Since  peak  reduced  electric 
field  in  these  calculations  varied  from  (T/7V)peak  ~  600  Td  (at  p  =  30  Torr)  to  (E/N)pea k  ~  200  Td  (at 
p  =  760  Torr),  vibrational  excitation  by  electron  impact  was  neglected.  For  E/N  >  200  Td,  discharge 
energy  fraction  loaded  into  the  vibrational  energy  mode  of  nitrogen  does  not  exceed  10%.  Heat  transfer 
from  the  gas  heated  in  the  discharge  was  also  neglected.  Using  sensitivity  analysis,  a  reduced  kinetic  model 
incorporating  15  species  (N2,  02,  O,  03,  NO,  N,  0(XD),  N2(A3F),  N2(B3n),  N2(a/:L£),  N2(C3n),  e-,  N+, 
0+,  Oj )  and  42  processes  was  obtained  from  the  full  air  plasma  model  including  46  species  and  395  processes 
(see  Tables  1-2).  Since,  at  the  high  peak  reduced-electric- fields  involved  here,  the  rate  of  electron  impact 
ionization  greatly  exceeds  the  rate  of  electron  attachment,  processes  involving  negative  ions  do  not  affect  the 
energy  balance.  During  the  sensitivity  analysis,  the  main  criterion  used  was  the  effect  of  individual  processes 
on  time-dependent  energy  fraction  thermalized  after  the  discharge  pulse.  The  results  of  calculations  for 
p  =  30  Torr  are  summarized  in  Figs.  3-6,  and  for  p  =  760  Torr  in  Figs.  7-10. 

From  Fig.  3,  it  can  be  seen  that  at  p  =  30  Torr,  peak  reduced  electric  field  and  peak  ionization  fraction 
are  approximately  (E/N)pea k  =  580  Td  and  (ne/7V)pea k  =  8  x  10-4  {ne  =  7.7  x  1014  cm-3),  respectively. 
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Figure  4  shows  that  energy  fraction  thermalized  during  the  first  10  ps  after  the  pulse  is  approximately 
31%.  Comparing  Figs.  4  and  5  demonstrates  that  the  energy  thermalization  on  this  time  scale  is  primarily 
due  to  collisional  quenching  of  N2(A3£),  N2(B3II),  N2(a/:L£),  N2(C3n),  and  0(1D),  as  well  as  electron-ion 
recombination.  Further  energy  thermalization,  on  a  much  longer  time  scale  of  10  ps  to  100  ms,  occurs 
primarily  in  reactions  of  nitrogen  and  oxygen  atoms  with  oxygen  molecules,  N  +  02^N0  +  0  and  O  + 
02  +  M  -A  O3  +  M,  resulting  in  nitric  oxide  and  ozone  formation,  as  well  as  in  reaction  O  +  O3  — >  02  + 
02.  Figure  4  also  demonstrates  the  accuracy  of  the  reduced  kinetic  model  as  compared  to  the  full  model. 

The  results  obtained  at  p  =  760  Torr  are  shown  in  Figs.  7-10.  In  this  case,  both  peak  reduced  electric 
field  and  peak  ionization  fraction  are  lower,  (E/lV)peak  =  200  Td  and  (ne/7V)peak  =  2  x  10-4  (ne  = 
4.9  x  1015  cm-3),  respectively  (see  Fig.  7).  From  Fig.  8,  it  can  be  seen  that  energy  fraction  thermalized 
within  1  ps  after  the  pulse  is  approximately  23%.  Again,  energy  thermalization  on  this  time  scale  is  due 
to  collisional  quenching  of  N2(A3E),  N2(B3II),  N2(a/:LI1),  N2(C3II),  and  0(1D)  (compare  Figs.  8  and  9), 
followed  by  an  additional  25%  of  input  energy  thermalized  on  a  long  time  scale,  1-300  ps  after  the  pulse,  in 
chemical  reactions  of  nitric  oxide  and  ozone  formation.  The  reduced  kinetic  model  obtained  based  on  the 
modeling  results  at  p  =  30  Torr  remains  accurate  (see  Fig.  8). 

The  rest  of  the  zero-dimensional  modeling  calculations  were  consistent  with  the  low-pressure  and  the 
high-pressure  cases  shown  in  Figs.  3-10.  For  the  entire  pressure  range  tested,  p  =  30  —  760  Torr,  the 
rapidly  thermalized  discharge  energy  fraction  was  in  the  range  of  25-30%.  The  time  scale  for  rapid  energy 
thermalization,  predicted  by  the  model,  varies  from  2-3  ps  at  p  =  30  Torr  to  0.2-0. 3  ps  at  p  =  760  Torr  (e.g., 
see  Figs.  4  and  8).  Thus,  rapid  heating  after  a  discharge  pulse  coupling  100  meV/molecule  to  dry  air  would 
result  in  localized  heating  by  80-100  K. 

Rate  coefficients  of  electron  impact  processes  incorporated  into  the  reduced  kinetic  model,  predicted  by 
the  Boltzmann  equation  solver  and  approximated  as  functions  of  the  reduced  electric  field,  are  listed  in 
Table  2.  These  rates,  as  well  as  rate  coefficients  of  chemical  reactions,  excited  electronic  species  quenching, 
and  electron-ion  recombination,  have  been  used  in  the  coupled  electric  discharge  /  compressible  flow  model 
to  predict  heating  and  compression  wave  formation  in  a  nanosecond-pulse  discharge  in  a  geometry  with 
planar  symmetry. 


III.  Discharge  Computations 

The  15-species  reduced  kinetic  model  was  applied  to  transient,  one-dimensional  discharge  computations. 
The  following  models  were  employed  to  predict  particle  motion:  a  drift-diffusion  formulation  for  the  charged 
particles,  a  diffusion  equation  for  the  neutrals,  and  a  mass- averaged  fluid  formulation  for  the  bulk  gas.  The 
Poisson  equation  was  solved  for  the  electric  potential.  A  detailed  description  of  the  overall  physical  model 
is  given  below. 

For  computational  efficiency,  the  calculations  were  carried  out  in  two  stages.  The  first  stage  encompassed 
the  period  in  which  an  external  waveform  was  applied  to  the  electrodes.  For  this  stage,  the  full  physical 
model  was  employed.  Shortly  after  the  external  potential  was  turned  off,  electromagnetic  effects  and  charged 
particle  motion  became  negligible.  In  this  second  stage  of  the  problem,  we  set  the  electric  field  to  zero,  and 
imposed  neutrality  by  setting  the  electron  number  density  to  an  appropriate  value. 

A.  Governing  Equations 

A  drift-diffusion  formulation  was  used  for  the  charged  particles: 

din 

+V-(nsw  ±ns/isE  -  DsVns)  =  los  (1) 

and  a  diffusion  formulation  was  used  for  the  neutrals: 

dm 

-7^  +  V-  (nsw  -  DsX7ns)  =  cus  (2) 

Here  the  electric  field  is  E  =  — V0,  where  0  is  the  electric  potential.  The  subscript  s  indicates  the  species 
number.  For  each  species-8,  the  number  density  is  ns,  the  mobility  is  ps,  the  diffusion  coefficient  is  Ds ,  and 
the  rate  of  production  in  reactions  is  ujs.  The  bulk  fluid  velocity  is  w. 
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A  mass-averaged  fluid  dynamic  formulation  was  used  to  model  the  motion  of  the  gas  as  a  whole: 


^+V-(pw)=0 

d 

—  (pw)  +  V-(/oww  +pl)  =  V-t  +  CE  (3) 

at 

d 

—  [/o(e  +  w2/2)}  +  V- [pw(e  +  w2/2 )  +pw]  =  V-(r-w  -  Q)  +  Ej 

Here  p  is  the  mass  density,  p  is  the  pressure,  r  is  the  viscous  stress  tensor,  £  is  the  space  charge,  j  is  the 
total  electric  current,  e  is  the  internal  energy,  Q  is  the  heat  flux.  It  useful  to  also  define  the  dissipated  power 
P  =  E-J.  where  J  =  j  —  (w. 

The  relationship  between  the  species  properties  and  the  overall  properties  are  as  follows: 


P  =  yprnsns  p  = 

S 

pe  =  y  msnses  pw  = 

S 

C  =  Y1  <isns  i  = 

s 

where  the  mass  per  particle  is  ras,  the  charge  per  particle 
the  species  internal  energy  is  es  =  h®  +  CV^ST.  Further, 
w  ±  psE  -  (DsWns)/ns. 

The  viscous  stress  and  heat  flux  are  as  follows: 

T  =  pv  [(Vw)  +  (Vw)T  -  §V-wl] 

Q  =  -fcVT  +  psUshs  (5) 

where  pv  is  the  viscosity,  k  is  thermal  conductivity,  Us  =  vs  —  w  is  the  diffusion  velocity,  and  hs  =  es-\-ps/ps 
is  the  species  enthalpy. 

Given  the  very  high  diffusion  velocities  in  the  cathode  layer  of  a  nanosecond-pulse  discharge,  careful 
consideration  must  be  made  of  the  validity  of  a  formulation  based  on  the  mass-averaged  global  conservation 
laws.  The  mass  averaged  formulation  assumes,  for  example,  that  the  species  internal  energy  is  the  same  for 
the  mass-averaged  reference  frame  as  for  the  species  reference  frame.  In  the  Appendix,  we  show  that  the 
formulation  used  here  is  valid  as  long  as  the  kinetic  energy  associated  with  diffusion  is  small  compared  to 
the  mixture  internal  energy.  For  the  present  work,  this  ratio  is  at  most  ^nsmsUg /(pe)  ~  10— 3 ,  so  the  model 
is  well  justified. 

The  Poisson  equation  is  employed  to  compute  the  electric  potential: 

V20  =  -C/eo  (6) 


J2Ps 

S 

msns\rs  ^ 

S 

Y,  qsnsvs 

S 

is  gs,  the  partial  pressure  is  ps  =  nskBT ,  and 
the  laboratory-frame  species  velocity  is  vs  = 


where  eo  is  the  permittivity  of  free  space. 

B.  Gas  Properties  and  Boundary  Conditions 

The  15-species,  42-process  formulation  described  previously  was  employed,  including  neutrals,  ions,  electrons, 
and  electronically-excited  species  (Tables  1-2).  Data  for  mobility  and  diffusion  coefficient  for  each  species 
of  heavy  particle  were  take  from  the  literature. 13,41-43  Correlations  of  electron  temperature  and  electron 
mobility  with  reduced  electric  field  (Fig.  11)  were  developed  from  the  Boltzmann  equation  solutions  described 
in  Sec.  II.  The  viscosity  and  thermal  conductivity  of  the  bulk  gas  were  based  on  standard  correlations  for 

•  44 

air. 

A  one-dimensional  computational  domain  was  employed  to  represent  a  single  dielectric  barrier  discharge 
(Fig.  12).  No-slip  boundary  conditions  with  a  constant  temperature  wall  were  employed  for  the  bulk  gas.  A 
zero  wall- normal  derivative  was  imposed  for  the  neutral  species. 

Standard  boundary  conditions  were  employed  for  the  charged  particles.  First,  the  conditions  at  the  wall 
were  determined  by  setting  the  normal  derivative  to  zero.  Then,  if  the  provisional  ion  flow  was  away  from 
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the  boundary,  the  ion  flux  was  set  to  zero.  The  electron  number  density  was  determined  through  a  balance 
between  the  kinetic  flux  to  the  wall  and  secondary  emission.  For  certain  cases,  these  boundary  conditions 
led  to  numerical  instability  at  the  exposed  electrode  boundary  (anode).  When  this  occurred,  the  simplified 
boundary  conditions  described  in  Ref.  45  were  employed. 

For  the  bare  electrode  the  potential  was  specified  as  zero  (grounded).  An  alternative  boundary  condition 
was  employed  for  the  powered  electrode  with  a  dielectric  coating  (see  also  Ref.  46).  The  dielectric  layer 
was  assumed  to  be  sufficiently  thin  that  a  linear  potential  profile  (uniform  electric  field  E^)  was  a  good 
approximation.  The  electric  field  inside  the  dielectric  was  related  to  the  electric  field  E  at  the  surface  through 
the  relation  eoE  — ereoE^  =  an,  where  a  is  the  surface  charge  density  and  n  is  the  unit  normal  vector  pointing 
into  the  computational  domain.  The  surface  charge  was  determined  by  integrating  dcr/dt  =  —  j  •  n  for  each 
surface  point,  using  a  time- marching  scheme  analogous  to  that  of  the  main  governing  equations. 

C.  Numerical  Methods 

The  calculations  were  carried  out  using  the  Air  Force  Research  Laboratory  code  HOPS  (Higher  Order  Plasma 
Solver).47  50  The  code  includes  several  physical  models  and  numerical  schemes.  Here,  the  physical  model 
consisting  of  Eqs.  (l)-(3)  and  (6)  was  solved  using  an  implicit,  second-order,  upwind  formulation.  All  the 
equations  were  solved  in  a  nondimensional  form  that  has  been  described  in  previous  papers.47  50 

Time  integration  of  the  conservation  equations  (l)-(3)  was  carried  out  using  a  second-order  implicit 
scheme,  based  on  a  three-point  backward  difference  of  the  time  terms.  The  formulation  is  similar  to  the 
standard  technique  of  Beam  and  Warming,51  but  adapted  to  a  multi-fluid  formulation  with  different  models 
for  particle  motion. 

Approximate  factoring  and  quasi-Newton  subiterations  were  employed.  The  implicit  terms  were  linearized 
in  the  standard  Thin  layer’  manner.  The  implicit  terms  were  evaluated  with  second-order  spatial  accuracy, 
yielding  a  block  tridiagonal  system  of  equations  for  each  factor.  The  species  were  loosely  coupled,  limiting  the 
rank  of  the  flux  Jacobian  matrices  to  the  order  of  the  moment  model  (one  for  the  drift-diffusion  formulation, 
five  for  the  overall  conservation  equations).  Each  factor  was  solved  in  turn  using  a  standard  block  tridiagonal 
solver,  and  the  change  in  the  solution  vector  of  conserved  variables  was  driven  to  zero  by  the  subiteration 
procedure  at  each  time  step.  Three  applications  of  the  flow  solver  per  time-step  were  employed  for  the 
present  work. 

For  the  mass-averaged  fluid  model  employed  for  the  bulk  gas,  the  Roe  scheme52,53  was  employed  for 
the  inviscid  fluxes.  For  the  drift-diffusion  model,  a  simple  upwinding  scheme  was  employed,  based  on  the 
convection-drift  velocity.  This  is  similar  to  the  approach  of  Surzhikov  and  Shang.45  In  both  formulations, 
stability  was  enforced  using  the  minmod  limiter  in  the  MUSCL  formalism.54 

The  Poisson  equation  (6)  was  solved  at  the  end  of  each  sub-iteration  in  the  implicit  time-marching  scheme. 
(Using  this  strategy  with  the  subiteration  procedure  gives  about  the  same  improvement  in  stable  time-step 
as  methods  based  on  the  linearization  of  the  right-hand-side  of  the  Poisson  equation.55)  An  approximately 
factored  implicit  scheme  was  employed,  adapted  from  the  approach  described  by  Holst.56  The  formulation 
of  the  implicit  scheme  was  analogous  to  that  of  the  conservation  equations,  with  linearization  of  the  implicit 
terms,  approximate  factoring,  and  an  iterative  procedure  that  drives  the  change  in  the  solution  to  zero.  The 
spatial  derivatives  were  evaluated  using  second-order  central  differences,  and  the  system  was  solved  using 
the  Thomas  tridiagonal  algorithm.57 

As  mentioned  earlier,  the  calculations  were  carried  out  in  two  phases.  The  first  stage  of  the  calculations 
encompassed  the  first  100  ns  of  the  discharge,  and  employed  the  full  physical  model  discussed  above,  Eqs.  (1)- 
(3)  and  (6).  Since  electromagnetic  effects  and  charged  particle  motion  became  negligible  after  the  input  pulse 
died  away  (after  about  24  ns),  in  the  second  stage  of  the  computations  (0.1  ps  to  100.1  ps),  the  electric 
field  was  set  to  zero  and  neutrality  was  enforced  by  appropriately  setting  the  electron  number  density.  This 
approach  resulted  in  a  substantial  savings  in  computational  cost. 

For  the  calculations  presented  here,  a  uniform  grid  of  1001  points  across  the  gap  was  employed.  Grid 
resolution  studies  presented  in  a  previous  paper28  indicate  that  this  level  of  grid  resolution  is  sufficient  for 
this  problem.  The  time  step  used  for  the  Phase  1  calculations  was  1  ps,  and  the  time-step  for  Phase  2  was 
5  ns. 
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D.  Results 


In  the  experiments  illustrated  in  Figs.  1-2,  typical  freestream  conditions  for  the  neutral  gas  were  a  speed 
of  715  m/s,  a  pressure  of  159.5  Pa,  and  a  temperature  of  56  K.  Sample  calculations  are  presented  here 
for  a  one-dimensional  discharge  under  conditions  representative  of  the  stagnation  region  of  the  cylinder 
flow  experiment.  The  corresponding  stagnation  conditions,  computed  using  the  Rayleigh  supersonic  Pitot 
formula58  were  4.74  kPa  (36  Torr)  and  310.3  K.  For  each  case  in  the  calculations  discussed  below,  the  initial, 
uniform  state  of  the  neutral  gas  was  set  to  these  stagnation  values. 

The  configuration  considered  here  is  illustrated  in  Fig.  12.  The  problem  is  one-dimensional.  In  the 
simulations,  the  right  electrode  was  grounded  and  the  left  electrode,  which  was  coated  with  a  thin  dielectric 
layer,  was  powered  with  the  input  signal  Vs  =  —Vo  exp  [—  (t  —  to)2/r2] .  For  the  present  calculations,  we 
assumed  Vo  =  20  kV,  r  =  3  ns,  and  to  =  12  ns. 

The  discharge  gap  was  taken  to  be  10  mm.  The  initial  mole  fraction  of  the  electrons  and  each  of  the 
neutral  minor  species  was  taken  to  be  1  x  10-10.  The  mole  fraction  for  each  ion  species  was  equal,  and  set 
so  that  the  space  charge  was  zero.  The  initial  electric  field  was  zero.  The  dielectric  coating  was  assumed  to 
be  2  mm  thick,  with  a  relative  dielectric  constant  of  er  =  3.8,  chosen  to  be  representative  of  fused  quartz. 
The  secondary  emission  coefficient  was  7sem  =  0.05. 

A  previous  paper13  presented  an  analytical  model  of  nanosecond-pulse  dielectric  barrier  discharges,  in 
which  an  electron-free  sheath  region  was  coupled  to  a  uniform,  quasi-neutral  plasma  region.  The  model 
provided  closed- form  expressions  for  the  discharge  properties,  and  displayed  good  agreement  with  experiment. 
We  employ  this  model  here  to  provide  a  basic  check  of  the  correctness  of  our  computer  code,  and  to  introduce 
the  general  features  of  the  nanosecond-pulse  discharge  in  one  dimension.  Figure  13  compares  the  discharge 
properties  in  the  plasma  predicted  by  this  analytical  model  and  by  the  one-dimensional  computations.  Note 
that  to  =  100  ns  for  the  analytical  model,  but  to  =  12  ns  for  the  numerical  computations. 

Figures  13a-b  show  the  applied  electric  field  (field  in  the  absence  of  plasma),  the  computed  field,  and 
the  ionization  fraction.  Some  quantitative  differences  are  present,  but  the  general  shape  of  the  curves  is 
in  agreement.  Because  of  the  shielding  effect  of  space  charge  in  the  cathode  sheath  and  surface  charge 
accumulation  on  the  dielectric,  the  field  in  the  plasma  is  less  than  the  applied  field.  The  one-dimensional 
computations  predict  a  peak  electric  held  at  9.4  ns,  before  the  input  pulse  maximum  at  12  ns,  followed  by 
a  smaller  peak  of  opposite  sign  at  14.5  ns.  The  peak  electron  mole  fraction  reaches  about  2  x  10-6. 

Figures  13c-d  show  the  reduced  electric  held  E/N ,  the  dissipated  power  EJ,  and  the  increase  in  gas 
energy  density  A (pe)  (including  chemical  and  thermal  energy).  Again,  the  results  are  qualitatively  consistent 
between  the  two  models,  with  moderate  quantitative  differences. 

For  the  one-dimensional  numerical  computations,  the  peak  reduced  electric  held  is  about  750  Td.  Peak 
dissipated  power  (1.0  x  1011  W/m3)  occurs  as  the  reduced  electric  held  in  the  plasma  is  rapidly  falling,  but 
before  the  peak  in  the  applied  electric  held.  The  analytical  model  predicts  a  total  energy  transferred  to 
the  gas  of  about  0.70  meV/molecule,  whereas  the  numerical  computations  predict  about  one-third  of  this 
value  (0.2  meV/molecule  or  38  J/m3).  This  difference  is  primarily  due  to  the  use  of  different  expressions  for 
the  ionization  rate  coefficient  as  a  function  of  reduced  electric  held  in  the  analytical  model  and  the  present 
numerical  computations.  If  the  same  expressions  for  ionization  rate  are  used  in  the  analytical  model  and 
the  numerical  simulations,  the  agreement  for  the  predicted  pulse  energy  coupled  to  the  plasma  improves 
considerably  (see  Ref.  13). 

Recall  that  for  the  corresponding  zero-dimensional  computations  discussed  in  Sec.  II  (Figs.  3-6),  a  dis¬ 
charge  energy  loading  on  the  order  of  100  meV/molecule  was  required  to  replicate  the  experimental  results, 
with  about  30  meV/molecule  thermalized  with  0.1  ms.  In  the  present  one-dimensional  calculations,  we 
obtained  a  discharge  energy  loading  of  about  0.2  meV/molecule  in  the  plasma.  The  plasma  electric  held  is 
comparable  for  the  two  the  formulations  (580  Td  for  zero-dimensional  vs.  750  Td  for  one-dimensional),  but 
the  ionization  fraction  in  the  plasma  is  two  orders  of  magnitude  higher  for  the  zero-dimensional  calculations 
(8  x  1(T4  vs.  2  x  1(T6). 

The  discrepancy  between  the  predictions  of  the  two  models  occurs  because  of  plasma  self-shielding  in 
the  one-dimensional  calculations.  This  effect  rapidly  reduces  the  electric  held  in  the  plasma  and  limits  the 
peak  electron  density.  The  results  demonstrate  that  using  zero-dimensional  plasma  kinetic  modeling,  with 
an  imposed  voltage  waveform,  is  not  appropriate  even  in  a  simple,  symmetric  configuration,  since  it  greatly 
overestimates  both  peak  electron  density  and  energy  coupled  to  the  plasma. 

Figure  14  shows  discharge  prohles  at  several  stages  in  the  computations.  The  species  are  grouped: 
the  prohles  of  the  total  number  densities  of  all  ions  (Nj  ,  O^,  and  0+),  electrons  (e-),  excited  neutrals 
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(N2(A3E),  N2(B3n),  N2(a/1S),  N2(C3II),  and  0(1D)),  and  ground-state  neutrals  (0,  O3,  NO,  N)  are  shown. 
The  concentrations  of  N2  and  02  remain  almost  constant,  and  are  omitted.  The  onset  of  breakdown  and 
appearance  of  the  peak  reduced  electric  field  are  illustrated  in  Figs.  14a-b  (9-10  ns).  The  sheath  forms 
rapidly  at  the  left  (covered)  electrode,  is  about  1  mm  thick,  and  carries  most  of  the  potential  drop  across  the 
gap.  Most  of  the  production  of  new  species  occurs  near  the  sheath  edge  and  in  the  plasma.  At  later  stages 
in  the  computation  (Fig.  14c,  14  ns),  the  electric  field  reverses  sign  in  the  plasma  (see  also  Ref.  21),  with  an 
accompanying  local  maximum  in  time  of  the  field  magnitude.  At  the  end  of  first  state  of  the  calculations 
(Fig.  14d,  100  ns),  quasi- neutrality  prevails  in  the  domain,  and  the  electric  field  is  negligible.  Since  diffusion 
has  not  had  time  to  act,  a  number  of  local  maxima  are  present  in  the  species  number  density  profiles.  These 
appear  to  reflect  the  motions  of  the  sheath  edge. 

Figure  15  shows  the  time-history  in  the  plasma  (x  =  8  mm)  of  the  mole  fraction  of  each  of  the  species 
except  N2  and  02,  which  remain  essentially  constant.  Results  are  shown  for  both  the  Phase  1  computations 
(0-0.1  ps,  space  charge  included)  and  Phase  2  computations  (0.1  ps-100.1  ps,  zero  electric  field  and  neutrality 
imposed).  The  concentrations  of  charged  and  excited  particles  increase  suddenly  when  the  electric  field 
becomes  strong,  and  decrease  gradually  after  the  input  waveform  ends.  The  longest-lived  species  are  Oj 
and  N2(A3E).  The  radicals  N  and  O  appear  at  breakdown,  whereas  NO  and  O3  increase  gradually  after  the 
end  of  the  input  waveform. 

We  see  from  Fig.  15  that  breakdown  results  in  the  rapid  formation  of  new  species,  which  gradually 
recombine  over  a  longer  time  scale.  Figure  16  addresses  the  accompanying  energy  redistribution.  For  a 
station  in  the  quasi-neutral  plasma  (x  =  8  mm),  the  distribution  of  chemical  energy  over  different  groups  of 
species  is  shown  in  Fig.  16a,  and  the  distribution  between  chemical  and  thermal  energy  is  shown  in  Fig.  16b. 
Corresponding  plots  for  a  station  near  the  sheath  edge  (x  =  1  mm),  are  shown  in  Figs.  16c-d. 

Through  the  dissipative  power  term  E  J,  breakdown  converts  part  of  the  input  electric  energy  into 
chemical  energy  of  new  species  over  a  time  scale  of  about  1  ns,  then  recombination  reactions  convert  the 
stored  chemical  energy  to  thermal  energy  over  roughly  a  10  ps  time  scale.  Initially,  most  of  the  chemical 
energy  is  stored  in  charged  particles  and  excited  neutrals.  With  time,  this  energy  thermalizes  through 
electron-ion  recombination  and  collisional  quenching  of  electronically  excited  species.  At  the  same  time 
additional  energy  goes  into  neutral  species  in  the  ground  electronic  states  (N,  O,  NO,  and  O3),  tending  to 
reduce  the  efficiency  of  the  device  in  rapidly  converting  electrical  energy  to  heat. 

In  the  plasma,  most  of  the  input  energy  is  initially  stored  in  chemical  energy,  whereas  near  the  sheath 
edge  this  fraction  is  relatively  small.  In  ongoing  work,  we  are  investigating  how  nonlocal  effects  and  the 
uncertainty  in  reaction  rates  in  this  region  of  very  high  reduced  electric  field  influence  this  result. 

The  oscillations  in  internal  energy  seen  in  Fig.  16c  for  times  of  10-20  ns  are  related  to  the  rapid  motion 
of  the  sheath  edge  during  that  time  interval.  (See  Fig.  14.) 

Also  shown  in  Figs.  16b  and  16d  is  a  comparison  of  the  time-integral  of  the  input  power  density  E-  J  dt 
to  the  change  in  total  internal  energy.  For  the  plasma  region,  these  are  seen  to  agree  almost  exactly,  indicating 
that  the  heating  is  essentially  a  zero-dimensional  phenomenon:  spatial  gradient  terms  like  convection  and 
heat  conduction  act  over  much  longer  time  scales.  For  the  sheath  region,  however,  fluid  motion  is  seen  to 
begin  to  carry  away  thermal  energy  within  100  ns. 

The  acoustic  time  scale  mentioned  previously  is  nicely  illustrated  in  Figs.  16b  and  16d:  after  several 
microseconds,  acoustic  waves  form  due  to  the  rapid  heating,  creating  the  spikes  apparent  in  the  plots.  Only 
energy  that  is  released  as  heat  on  a  time  scale  shorter  than  the  acoustic  time  can  contribute  to  the  formation 
of  these  waves. 

The  acoustic  waves  are  more  apparent  in  Fig.  17,  which  shows  the  heating  and  gas  motion  induced  by 
the  discharge.  The  majority  of  the  heating  occurs  near  the  edge  of  the  cathode  sheath,  forming  weak  waves 
that  travel  across  the  domain.  For  example,  consider  the  velocity  profiles  in  Fig.  17b,  and  note  the  wave 
structure.  Two  waves,  traveling  in  opposite  directions,  appear  near  x  —  l  mm  for  the  0.1  ps  profile  (red 
curve).  After  4.1  ps  (green  curve),  the  left-running  wave  has  reflected  off  the  left  boundary  and  trails  the 
other  compression  wave.  Note  also  the  growth  in  wave  strength  between  these  two  times;  the  waves  are 
driven  by  exothermic  reactions  in  a  manner  analogous  to  a  detonation. 

For  this  case,  the  waves  are  relatively  weak,  with  a  peak  gas  velocity  of  about  3  m/s  and  peak  temperature 
rise  of  about  11  K.  The  stationary  background  gas  experiences  less  than  a  1  K  temperature  rise  over  the 
time  covered  in  the  simulation.  For  an  air  temperature  of  310.8  K,  the  speed  of  sound  is  about  353  m/s. 
Estimating  the  wave  speed  by  measuring  displacements  in  Fig.  17a,  we  find  355  m/s,  which  is  very  slightly 
supersonic. 
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Substantially  supersonic  wave  speeds  (up  to  480  m/s)  have  been  observed  experimentally.20  In  ongo¬ 
ing  work,  we  are  attempting  to  replicate  this  stronger  gasdynamic  interaction  by  altering  the  simulation 
conditions  to  increase  the  coupled  energy  density. 

IV.  Summary  and  Conclusions 

Numerical  calculations  were  carried  out  to  examine  the  physics  of  the  operation  of  a  nanosecond-pulse, 
single  dielectric  barrier  discharge  in  a  configuration  with  planar  symmetry.  This  simplified  configuration  was 
chosen  as  a  vehicle  to  develop  a  physics-based  nanosecond  discharge  model,  including  realistic  air  plasma 
chemistry  and  compressible  bulk  gas  flow.  Discharge  parameters  (temperature,  pressure,  and  input  wave¬ 
form)  were  selected  to  be  representative  of  recent  experiments  on  bow  shock  control  with  a  nanosecond 
discharge  in  a  Mach  5  cylinder  flow. 

First,  a  reduced  plasma  kinetic  model  (15  species  and  42  processes)  was  developed  by  carrying  out  a 
sensitivity  analysis  of  zero-dimensional  plasma  computations  with  an  extended  chemical  kinetic  model  (46 
species  and  395  processes).  Transient,  one-dimensional  discharge  computations  were  then  carried  out  using 
the  reduced  kinetic  model,  incorporating  a  drift-diffusion  formulation  for  each  species,  a  self-consistent  com¬ 
putation  of  the  electric  potential  using  the  Poisson  equation,  and  a  mass-averaged  gas  dynamic  formulation 
for  the  bulk  gas  motion. 

A  grid  converged  solution  and  reasonable  comparison  to  a  validated  analytical  model  indicate  that  the 
computations  reflect  an  accurate  solution  of  the  mathematical  model.  The  computational  results  qualita¬ 
tively  reproduce  many  of  the  features  observed  in  experiments,  including  the  rapid  thermalization  of  the 
input  electrical  energy  and  the  consequent  formation  of  a  weak  shock  wave.  The  results  illustrate  how  input 
electrical  energy  is  rapidly  transformed  (over  roughly  1  ns)  at  breakdown  into  ionization  products,  dissoci¬ 
ation  products,  and  electronically  excited  particles,  and  how  thermalization  occurs  over  a  relatively  longer 
time-scale  (roughly  10  ps). 

The  motivation  for  the  present  work  is  modeling  nanosecond-pulse,  dielectric  barrier  discharges  for  ap¬ 
plications  in  high-speed  flow  control.  The  effectiveness  of  such  devices  as  flow  control  actuators  depends 
crucially  on  the  rapid  thermalization  of  the  input  electrical  energy,  and  in  particular  on  the  rate  of  quench¬ 
ing  of  excited  electronic  states  of  nitrogen  molecules  and  oxygen  atoms  and  on  the  rate  of  electron-ion 
recombination. 

Future  work  will  include  using  the  coupled  nanosecond  discharge  /  compressible  flow  model  developed 
in  the  present  work  for  simulation  of  surface  nanosecond  pulse  discharges.  The  main  difference  of  the 
present  approach  from  other  recent  studies  of  nanosecond  pulse  discharges  in  air24  is  that  the  present  model 
incorporates  the  kinetics  of  energy  storage  and  thermalization. 

In  addition,  various  flow  control  applications  will  be  explored.  One  option  under  consideration  is  replacing 
the  cylinder  model  with  a  blunt  oblique  shock  generator,  and  examining  whether  the  nanosecond  discharge 
can  affect  an  impinging  shock  /  boundary  layer  interaction.  There  is  a  strong  motivation  to  control  such 
interactions  for  supersonic  engine  inlet  applications.27 

Appendix:  Mass- Averaged  Conservation  Equations 

The  present  paper  employs  a  mass-averaged  formulation  of  the  conservation  laws  for  the  gas  as  a  whole, 
Eq.  (3).  Here  we  present  a  brief  derivation  of  the  mass-averaged  formulation  from  the  conservation  laws  for 
individual  species,  and  discuss  the  range  of  applicability  of  the  conventional  mass- averaged  formulation  for 
our  applications.  (See  also  Ref.  59.) 

A  derivation  of  the  conservation  laws  for  the  individual  species,  as  moments  of  the  Boltzmann  equation, 
was  presented  in  a  previous  paper,49  and  is  also  addressed  elsewhere.60  62  Briefly,  the  moment  equations  are 
obtained  by  multiplying  the  Boltzmann  equation  by  a  conserved  quantity,  and  integrating  over  all  velocity 
space.  Considering  a  quantity  0(u)  associated  with  each  particle  of  species-8,  an  average  value  is  defined  as 
(0)s(x,  t)  =  f  0(u)fs(x,u,t)d3u,  where  where  u  is  the  particle  velocity,  fs  is  the  distribution  function,  and 
the  integral  is  over  all  velocity  space.  The  resulting  the  mass,  momentum,  and  energy  conservation  equations 
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for  species- s  are: 


dps 

dt 


+  V-(/9svs)  =  Ss 


d_ 

dt 


(psvs)  +  V-(psvsvs+psI)  =  V-rs  +psg  +  Cs  (E  +  vs  xB)  +  As 
d 

~Qt  [ps  (es  +  5^)]  +  V-  [psvs  (es  +  ^v2)  +psvs]  = 
V-[rs-vs  -  Qs]  +psvs-g  +  Csvs-E  +  Ms 


(7) 


where  the  mass  density  has  been  defined  as  ps  =  nsms  and  the  charge  density  as  (s  =  nsqs.  The  mean 
velocity  is  defined  as  vs  =  ( u)s  and  the  peculiar  velocity  as  Vs  =  u  —  vs.  The  source  terms  Ss ,  As,  and 
Ms  represent  the  exchange  between  species  of  particle  identity,  momentum,  and  energy  in  collisions.  Also 
defined  are  the: 

Kinetic  pressure:  ps  =  \nsms(y^) s 
Thermal  energy:  es  =  einM  +  \  {V^)S 
Viscous  stress:  rs  =  -nsms[(VsVs)s  -  |(VS2)SI] 

Heat  flux.  —  Qint,s  T  2^s Vs  )s 

The  terms  eint?s  and  Qint,s  represent  internal  molecular  energy  and  the  transfer  of  such  energy  by  diffusion 
(e.g.,  see  Ref.  62). 

The  conservation  of  mass,  momentum,  and  total  energy  in  collisions  requires  that  the  sum  over  all 
species  of  each  of  the  corresponding  collision  source  terms  is  zero:  Ss  =  0,  As  =  0,  Ms  =  0. 

These  properties  can  be  used  to  relate  the  species  conservation  laws  (7)  to  the  conservation  laws  for  the  gas 
as  a  whole. 

The  overall  density  is  defined  as  p  =  ^2sps-  Consider  the  alternative  peculiar  velocity  W  =  u  —  w, 
where  w  =  psvs/ P  is  the  mass- weighted  average  velocity.  The  species  diffusion  velocity  is  defined  as 
Us  =  (W)8  =  vs  —  w.  (Note  that  W  =  Us  +  Vs.)  In  the  mass-average  reference  frame,  the  moment  terms, 
analogous  to  Eq.  (8),  are: 

Ps  =  3 Ps (hE  )s  =  Ps  T  ^psUs 

is  =  Cint.s  +  \{W2)S  =  es+  \U2 

Ts  =  -Ps[{WW)s  -  l(W2)sl]  =TS-  ps[VsUs  -  |f/2I]  (9) 

Qs  =  Qint.a  +  \ps(WW2)s 

=  Q,  +  \ps(hs  +  ±U?)I-T3]-Va 


where  hs  =  es  +ps/ps  is  the  species  enthalpy.  Introducing  (9)  into  (7),  the  following  conservation  equations 
are  obtained: 


dps 

dt 


+  \7-(psw  +  ps\Js)  =  Ss 


+  PsUs)  +  V'(ftww  +  wpsus  +  psusw  +  psl)  = 

V-f s  +  psg  +  CSE  +  (CsW  +  CsUs)  x  B  +  As 

d 

^  [ps  (es  +  \w2)  +  psUs-w] 

+V-  [/9sw  (es  +  ±w2)  +/)SUS'WW  +  ps\Js\w2  +psw]  = 

Ts-w  -  QS1  +  (ps W  +  psUs)-g  +  E-(Csw  +  CsUs)  +  Ms 


V- 


Define  the  following  overall  properties: 


(10) 


e  =  ^2psis/p 

S 

p  =  EA 


C  =  ECa 


j  =  Ecsus 

s 

q  =  E^ 

s 


(ii) 
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The  terms  in  (11)  are,  respectively,  the  net  thermal  energy,  space  charge,  conduction  current,  pressure, 
viscous  stress,  and  heat  flux.  The  total  electrical  current  is  j  =  £w  +  J.  Using  these  new  definitions, 
and  summing  (10)  over  all  species,  the  global  conservation  laws  for  mass,  momentum,  and  total  energy  are 
obtained: 

^  +  V-(/3w)  =0 

Iw  +  V-tpww+pI)  =  V-r  +  pg  +  CE+j  xB 

d 

—  [p(e  +  \w2)\  +V-  [pw(e  +  \w2)  +  pw]  = 

V •  [r  •  w  -  Q]  +  pw-g  +  Ej 

where  all  the  collision  source  terms  have  summed  to  zero.  This  is  a  slightly  generalized  form  of  Eq.  (3). 

Note  the  appearance  in  Eq.  (9)  of  the  diffusion  velocity  terms  that  arise  from  the  change  of  reference 
frame  from  that  of  the  species  velocity  to  that  of  the  mass-averaged  velocity.  If  terms  of  order  U ^  are 
neglected,  all  the  quantities  except  the  heat  flux  are  the  same  in  the  mass-averaged  reference  frame  as  in  the 
species  reference  frame.  For  the  heat  flux  we  have  Qs  ~  Qs  +  psXJshs,  so  that  the  total  heat  flux  becomes 
Q  =  +  psXJshs  in  this  approximation,  which  is  the  form  used  in  Eq.  (5). 

The  error  introduced  by  this  approximation  in  the  summations  of  Eq.  (11)  is  small  if  \psU J  <C  pe,  in 
other  words  for  either  small  diffusion  velocities  or  small  mass  fractions.  This  is  a  good  approximation  for 
a  weakly-ionized  gas,  because  the  neutral  particles  have  high  number  density  but  small  diffusion  velocities, 
whereas  the  charged  particles  have  high  diffusion  velocities  but  low  number  densities. 

A  check  of  representative  values  of  \psU^ / {pe)  at  a  time  (10  ns)  close  the  the  peak  in  electric  field 
magnitude  in  the  present  calculations  illustrates  this  point.  In  the  numerical  solution,  this  parameter  for 
both  ions  and  electrons  has  value  of  about  1  x  10-6  in  the  plasma.  Values  for  the  electrons  drop  off  in  the 
sheath,  but  the  values  for  the  ions  peak  near  the  sheath  edge,  with  a  maximum  of  about  2  x  10-3  for  N^. 
Thus  the  mass- averaged  formulation  should  be  accurate  for  the  present  calculations,  but  a  careful  check  is 
warranted  if  calculations  are  carried  out  for  lower  density  discharges. 
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No. 

Species 

h°  (eV/particle) 

1 

n2 

0 

2 

02 

0 

3 

0 

2.58 

4 

03 

1.48 

5 

NO 

0.93 

6 

N 

4.90 

7 

O^D) 

2.0 

8 

N2(A3£) 

6.17 

9 

N2(B3n) 

7.35 

10 

3 

to 

8.40 

11 

N2(c3n) 

11.0 

12 

e_ 

0 

13 

n2+ 

16.36 

14 

0+ 

16.26 

15 

01 

12.70 

Table  1.  Species  incorporated  in  kinetic  model. 
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No. 

Reaction 

Rate  Expression 

Coefficients 

1 

N2  +  e_  — »  +  e_  +  e_ 

l°g10  k  =  A  +  B/(E/N) 

-7.76,  -37.0  (-7.17,  -61.1) 

2 

O2  +  c  — >  O^”  +  e  +  e 

log10  k  =  A  +  B/(E/N ) 

-8.34,  -30.7  (-7.42,  -71.0) 

3 

02  +  e“  ->  O  +  0+  +  e“  +  e“ 

log10  k  =  A  +  B/(E/N ) 

-7.94,  -32.2  (-7.78,  -35.2) 

4 

N2  +  e-  ->  N2(A3E)  +  e- 

log10fc  =  A+B/(E/N) 

-8.57,  -11.3  (-9.06,  15.3) 

5 

N2  +  e“  -»  N2(B3n)  +  e“ 

l°g10  k  =  A  +  B /{E/N) 

-7.97,  -13.2  (-8.14,  -3.42) 

6 

N2  +  e-  ->  N2(C3n)  +  e- 

log10  k  =  A  +  B/(E/N ) 

-7.82,  -22.0  (-7.45,  -37.0) 

7 

N2  +  e“  ->  N2(a,1E)  +  e“ 

log10  k  =  A  +  B/(E/N) 

-8.15,  -15.8  (-8.12,  -15.6) 

8 

N2  +  e“  ->►  N  +  N  +  e“ 

l°g10  k  =  A  +  B /{E/N) 

-7.96,  -22.7  (-7.34,  -50.2) 

9 

02  +  e“  ->  O  +  O  +  e“ 

l°g10  k  =  A  +  B /{E/N) 

-8.31,  -8.40  (-8.63,  7.10) 

10 

02  +  e-  ->  O  +  0(1D)  +  e- 

log10  k  =  A  +  B/(E/N ) 

-7.86,  -17.2  (-7.42,  -37.1) 

11 

N  +  02  ->■  NO  +  O 

k  =  ATn  exp {-Ea/T) 

0.110E-13,  1.00,  3150. 

12 

N  +  NO  -►  N2  +  O 

k  =  ATn 

0.105E-11,  0.50 

13 

O  +  O3  — >  O2  +  O2 

k  =  A  exp(— Ea/T) 

0.200E-10,  2300. 

14 

0  +  0  +  n2  ->  o2  +  n2 

k  =  A  exp(— Ea/T) 

0.276E-33,  -720. 

15 

0  +  0  +  o2  ->  o2  +  o2 

k  =  ATn 

0.245E-30,  -0.63 

16 

0  +  02  +  N2  — >  03  +  N2 

k  =  ATn 

0.558E-28,  -2.00 

17 

0  +  02  +  02  — >  03  +  02 

k  =  ATn 

0.861E-30,  -1.25 

18 

N2(A3E)  +  o2  ->•  n2  +  0  +  0 

k  =  0.170E-11 

19 

N2(A3E)  +  o2  -)•  n2  +  o2 

k  =  0.750E-12 

20 

N2(A3E)  +  O  -)•  N2  +  0(XD) 

k  =  0.300E-10 

21 

N2(A3E)  +  N2(A3S)  ->•  N2  +  N2(B3n) 

k  =  0.770E-10 

22 

n2(a3e)  +  n2(a3e)  ->  n2  +  N2(c3n) 

k  =  0.160E-09 

23 

N2(B3n)  +  N2  -»•  N2(A3E)  +  n2 

k  =  0.300E-10 

24 

N2(B3n)  ->  N2(A3E) 

k  =  0.150E+06 

25 

N2(B3n)  +  o2  n2  +  0  +  0 

k  =  0.300E-09 

26 

N2(a,1E)  +  N2  N2  +  N2 

k  =  0.200E-12 

27 

N2 (a,1E)  +  02  -4-  N2  +  O  +  0(XD) 

k  =  0.281E-10 

28 

N2(c3n)  +  n2  ->•  N2(B3n)  +  n2 

k  =  0.100E-10 

29 

N2(c3n)  -»•  N2(B3n) 

k  =  0.300E+08 

30 

N2(C3n)  +  02  ->  N2(A3E)  +  0  +  0 

k  =  0.301E-09 

31 

0(1D)  +  N2  -y  0  +  N2 

k  =  0.257E-10 

32 

0(1D)  +  02  -y  O  +  02 

k  =  0.400E-10 

33 

0+  +  02  -y  0+  +  O 

k  =  0.199E-10 

34 

N+  +  02  -y  N2  +  0+ 

k  =  0.600E-10 

35 

N+  +  e“  -y  N  +  N 

k  =  AT™ 

0.831E-05,  -0.50 

36 

Oj  +  e“  -y  O  +  O 

k  =  AT™ 

0.599E-04,  -1.00 

37 

Nj  +  e_  +  e_  — y  N2  +  e_ 

k  =  AT™ 

0.140E-07,  -4.50 

38 

O2"  +  c  +  c  — y  02  +  e 

k  =  AT™ 

0.140E-07,  -4.50 

39 

0+  +e_+e_^-0  +  e_ 

k  =  AT™ 

0.140E-07,  -4.50 

40 

N+  +  e"  +  M  -y  N2  +  M 

k  =  AT™ 

0.312E-22,  -1.50 

41 

Oj  +  e_  +  M  — y  02  +  M 

k  =  AT™ 

0.312E-22,  -1.50 

42 

0+  +e-  +  M— yO  +  M 

k  =  AT™ 

0.312E-22,  -1.50 

Table  2.  Reaction  mechanism.  Units  consistent  with  number  densities  in  cm-3,  and  temperatures  in  K,  reduced 
electric  field  in  1  x  10— 16  V  cm2  (10  Td),  one-body  rates  in  s— 1,  two-body  rates  in  cm3/s,  and  three-body  rates 
in  cm6/s.  Values  in  parentheses  are  alternative  curve  fits  for  E/N  >  50  X  10-16  V  cm2  ( E/N  >  500  Td). 
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(a)  Diagram  of  cylinder  mounted  in  (b)  Perspective  view  of  cylinder  model.  (c)  Side  view  of  cylinder  model, 

wind  tunnel. 


Figure  1.  Configuration  of  Mach  5  cylinder  flow  experiment. 


(a)  Time:  1  ps. 

Figure  2. 


(b)  Time:  2  ps.  (c)  Time:  3  ps.  (d)  Time:  4  ps. 

Side-view  schlieren  images  of  Mach  5  cylinder  flow  perturbed  by  pulsed  discharge. 
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Figure  3.  Reduced  electric  field  during  the  pulse;  ionization  fraction  predicted  by  the  full  and  reduced  air 
plasma  models.  Conditions:  p  =  30  Torr,  To  =  300  K,  discharge  energy  loading  100  meV/molecule. 


incV  molecule 


Figure  4.  Input  pulse  energy  and  energy  thermalized  after  the  pulse  predicted  by  the  full  and  reduced  air 
plasma  models  at  the  conditions  of  Fig.  3  (p  =  30  Torr,  Tq  =  300  K,  100  meV/molecule). 
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(a)  Full  kinetic  model. 


(b)  Reduced  kinetic  model. 


Figure  5.  Mole  fractions  of  dominant  neutral  species  at  the  conditions  of  Fig.  3  (p  =  30  Torr,  To  =  300  K 
100  meV/molecule). 


(a)  Full  kinetic  model. 


(b)  Reduced  kinetic  model. 


Figure  6.  Mole  fractions  of  dominant  charged  species  at  the  conditions  of  Fig.  3  (p  =  30  Torr,  To  =  300  K 
100  meV/molecule). 
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Figure  7.  Reduced  electric  field  during  the  pulse;  ionization  fraction  predicted  by  the  full  and  reduced  air 
plasma  models.  Conditions:  p  =  760  Torr,  To  =  300  K,  discharge  energy  loading  130  meV/molecule. 


meV  molecule 


Figure  8.  Input  pulse  energy  and  energy  thermalized  after  the  pulse  predicted  by  the  full  and  reduced  air 
plasma  models  at  the  conditions  of  Fig.  7  (p  =  760  Torr,  Tq  =  300  K,  130  meV/molecule). 
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(a)  Full  kinetic  model. 


(b)  Reduced  kinetic  model. 


Figure  9.  Mole  fractions  of  dominant  neutral  species  at  the  conditions  of  Fig.  7  (p  =  760  Torr,  To  =  300  K 
130  meV/molecule). 
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(a)  Full  kinetic  model. 


(b)  Reduced  kinetic  model. 


Figure  10.  Mole  fractions  of  dominant  charged  species  at  the  conditions  of  Fig.  7  (p  =  760  Torr,  To  =  300  K 
130  meV/molecule). 
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Electron  temperature.  eV 


Figure  11.  Electron  properties  as  a  function  of 
Boltzmann  equation  solutions. 


Electron  mobility,  cm2/Vs 


(b)  Electron  mobility  at  atmospheric  pressure  ( N  =  2.4  x 
1019  cm"3). 


electric  field.  Lines:  curve  fits;  points:  data  from 


Powered  Electrode 


Figure  12.  Diagram  of  computational  domain. 
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Electric  field.  V  ein  Ionization  fraction 


(a)  Electric  field  and  ionization  fraction,  analytical  model.  (b)  Electric  field  and  ionization  fraction,  1-D  computations 

{x  =  8  mm). 
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(c)  Reduced  electric  field,  power,  and  thermalized  energy,  an-  (d)  Reduced  electric  field  E/N ,  dissipated  power  E-  J,  and 
alytical  model.  change  in  internal  energy  A (pe),  1-D  computations  {x  = 

8  mm). 


Figure  13.  Comparison  of  analytical  model13  to  one-dimensional  computations  for  a  station  in  the  quasi-neutral 
plasma. 
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Figure  14.  Discharge  profiles  predicted  by  the  numerical  computations.  Ions:  N^"  ,  Oj,  and  0+;  electrons: 
e_;  excited  neutrals:  N2(A3£),  N2(B3II),  N2(a/1E),  N2(C3II),  and  0(1D);  ground-state  neutrals:  O,  O3,  NO, 
N. 
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Potential  (V)  Potential  (V) 


(b)  Electronically  excited  neutrals:  N2(A3£),  N2(B3I1), 

N2(a/:LE),  N2(C3n),  and  0(1D). 


(c)  Ground-state  neutrals:  O,  O3,  NO,  N. 


Figure  15.  Time  history  of  species  mole  fraction  in  the  plasma  (x  =  8  mm).  The  calculations  switch  from  the 
full  model  to  the  neutral,  zero-field  model  at  1  x  10_7s. 
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(a)  Distribution  of  chemical  energy  density  over  different 
groups  of  species.  (Plasma,  x  =  8  mm.) 


(c)  Distribution  of  chemical  energy  density  over  different 
groups  of  species.  (Sheath,  x  =  1  mm.) 


(b)  Distribution  of  energy  between  chemical  and  thermal 
modes.  (Plasma,  x  =  8  mm.)  Spikes  in  thermal  energy  corre¬ 
spond  to  acoustic  wave  motion. 


(d)  Distribution  of  energy  between  chemical  and  thermal 
modes.  (Sheath,  x  =  1  mm.)  Spikes  in  thermal  energy  corre¬ 
spond  to  acoustic  wave  motion. 


Figure  16.  Distribution  of  energy  in  the  plasma  and  in  the  sheath.  Charged  particles:  ,  Oj,  0+,  and  e  ; 

excited  neutrals:  N2(A3E),  N2(B3II),  N2(a/1E),  N2(C3II),  and  0(1D);  ground-state  neutrals:  O,  O3,  NO,  N; 
change  in  energy  density  A (pe),  integrated  power  f*  E-  J  di.  The  calculations  switch  from  the  full  model  to  the 
neutral,  zero-field  model  model  at  1  x  10— 7s. 
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Figure  17.  Profiles  of  the  properties  of  the  bulk  gas  for  selected  times  in  the  simulation.  Arrows  indicate  the 
direction  of  wave  motion. 
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